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We study two-dimensional soliton-soliton vector pairs in media with self-focusing nonlinearities 
and defocing cross-interactions. The general properties of the stationary states and their stability 
are investigated. The different scenarios of instability are observed using numerical simulations. 
The quasi-stable propagation regime of the high-power vector solitons is revealed. 
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I. INTRODUCTION 

Bright spatial solitons are localized in space self- 
induced structures that appear in various physical sys- 
tems as the result of a balance between nonlinear self- 
focusing and diffractive spreading As well known, in 
bulk Kerr media the two-dimensional (2D) spatial soli- 
tons are unstable and they either collapse or spread out 
depending on the power of the wave packet (see, e.g. 

The different mechanisms have been proposed to 
arrest the collapse, including higher-order linear Q and 
nonlinear JBl, @ effects, nonlocal response Q, orbital mo- 
mentum U and others [l[ , but in pure Kerr media for 
a single wave beam the problem of collapse is not yet 
solved. 

The multi-component self-trapped stationary struc- 
tures, known as vector solitons^ open the new possibil- 
ities for stabilization of over-critical beams. Vector soli- 
tons consist of more than one field components, which 
are supported not only by self-interaction of the same 
fields, but also by a cross-interactions between fields of 
different type. In nonlinear optical media the internal 
and intercomponent interactions are always of the same 
sign, so the localized soliton-soliton complexes either un- 
dergo the collapse [oj for focusing nonlinearities, or do 
not exist at all for defocusing nonlinearities. However, 
if media is self-attracting, but the interactions between 
different components become repulsive, the collapse ap- 
parently could be suppressed. Here the natural question 
arises of whether there is a realistic physical system with 
competing self-focusing and cross-defocusing nonlineari- 
ties? 

The first example of the media with appropriate non- 
linear properties is two-component BEG of ultracold 
atomic gases with Feschbach resonance management. 
The remarkable progress of experimental realization of 
BEC with tunable nonlinearities [lo| motivated inves- 
tigations of vector solitons with different signs of non- 
linearities [l^ - |22j . The soliton- vortex vector pair for 
two-component BEC with attractive intracomponent and 
repulsive intercomponent has been investigated in Ref. 
[23 ]. However, the ground state of this system is not 
studied yet. 

Second nonlinear physical system with tunable cross- 
interactions was found recently in plasma with bi-color 



laser beam Depending on a frequency difference, 

the cross-focusing or cross-defocusing is observed, while 
the self-interactions remain focusing. As was predicted in 
Ref. the balance between the competing nonlinear- 
ities can stabilize the system and result in a dynamical 
guiding of multi-color laser beam. 

It is remarkable that bi-color laser beam in plasma 
and matter-wave solitons in two-component BEC being 
two quite different physical systems, which belong to the 
opposite sides of the temperature scale, are described 
by the same model. This model is based on the set of 
two coupled nonlinear Schrodinger (NLS) equations with 
attractive internal and repulsive intercomponent cubic 
nonlinearities. In this paper we study the stationary so- 
lutions of the coupled NLS equations both numerically 
and analytically. Our variational treatment accounts for 
the essential modification of the soliton shape and agrees 
well with our numerical calculations. Stability of the ob- 
tained soliton-soliton pairs has been tested by numeri- 
cal simulations. We describe here the different scenarios 
of unstable evolution including aziumathally asymmetric 
modulational instability, which can substantially restrict 
propagation distance for high-power vector solitons. At 
the same time, we found out the condition for quasi- 
stable propagation of two-dimensional bright vector soli- 
tons. 



II. BASIC EQUATIONS 

Here we consider condition for the formation of self- 
induced structures and their stability on the basis of two 
coupled NLS equations: 

z^-l-(A^-M*i|2-l-a|vI/2|2)M/i=0, (1) 



dz 



i^ + (A^+a|*i|2 + |vI/2|2)vl/2 = 0, (2) 

where is a 2D Laplacian operator. The integrals of 
motion are the beam power (or number of particles for 
BEC) in each component 
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and the Hamiltonian 



where 




Also the momentum and the angular momentum are con- 
served, which we will not need in explicit form since these 
integrals vanish for all solutions under consideration. 

The model based on Eqs. (P), ^ is of broad physi- 
cal interest. In nonlinear optics these equations describe 
two noncoherently interacting wave beams propagating 
in z-direction, where z comes in units of Rayleigh length. 
The properties of bright vector solitons are well known for 
media with self-focusing and cross-focusing (cr > 0) non- 
linearities . While the parameter of coupling a varies 
over a broad range for different nonlinear optical media 
(depending on polarization state, nature of nonlinearity 
and anisotropy of media) , nevertheless the sign of cross- 
interaction coincides with the sign of self-interaction. 

The paraxial envelope equations ([T]), ^ have been de- 
rived in Ref. [ll| for the problem of bi-color laser beam 
propagation in plasmas. When a long bi-color laser beam 
with a frequency difference fl propagates in plasma, the 
relativistic cross-focusing provides the self-focusing ef- 
fect. At the same time, the cross-interaction is character- 
ized by the coupling constant a ~ (il^ — 2ujp)/{il'^ — w^), 
where cOp is the electron plasma frequency. The sign 
of the coupling parameter corresponds to cross-focusing 
(cr > 0, if < Wp, 17 > ^/2ujp), or cross-defocusing 
(cr < 0, if ojp < f2 < \/2ujp). The physical reason for 
such sharp tuning of the cross-interaction at fi ss Wp is 
that the ponderomotive force drives an electron plasma 
wave which acts as an either focusing or de-focusing chan- 
nel depending on the value of difference frequency. The 
analysis of the model given in Ref. [ll| was restricted 
by the semi-analytical method based on the dynamical 
equations for the parameters of bell-shaped Gaussian- 
type ansatzcs of both components. The suppression of 
the catastrophic self-focusing has been predicted in the 
frame of such a simplified treatment and supported by 
fully relativistic axially-symmetric PIC simulations of the 
laser beam dynamics over a large propagation distance. 

The NLS equations ([1]), ([2]), arc known also as Gross- 
Pitacvski equations, describe in mean-field approxima- 
tion the wave functions of two interacting BEG at ultra 
low temperature. The variable z should be replaced by 
dimensionless time t in context of BEG. Two equations 
correspond to two-component BEG of atoms of the same 
isotope in different hyperfine states. The atoms of BEG 
are trapped by strong planar external trap, and the order 
parameters in z-direction are frozen to the ground state, 
while the dynamics in (x, y) plane is described by Eqs. 
©, ©. 



III. STATIONARY SOLUTIONS 

We look for the radially-symmetric stationary solutions 
of Eqs. (H)) and ^ in the form 

vl/,(r,z)=^,(r)e'^^^ (5) 

where /3i and /32 are independent propagation constants. 
We are interested here in ground state solutions, so i'jir) 
can be treated as real functions. Let us choose as the 
scale of the radial coordinate and introduce the soliton 
parameter A = P2I Pi- Using the following scaling of the 
soliton profiles: ijji = y/]3iu{r), ■ip2 — \fWy-v(r) one can 
obtain the set of the stationary equations as follows 

— u -|- Arit -|- [x? -I- av^)u — 0, (6) 

- Xv + ArV + {(Tu'^ + v'^)v = 0, (7) 

where = + r ^® radial Laplacian. In this 
section we analyze both numerically and analytically two- 
parameter vector soliton families (with parameters A, and 
cr) for defocusing intercomponent nonlinearity (cr < 0). 

A. Numerical modelling 

The stationary equations have been solved by the sta- 
bilized iterative procedure similar to that described in 
Ref. The examples of radial profiles u{r) and v{r) 

for different coupling constants cr at fixed soliton param- 
eter A are given in Fig. [1] It is seen that for strong re- 
pulsive interaction between the two solitons their shapes 
change substantially. The field of first beam is squized 
out and forms the ring-like shell, while the second com- 
ponent is noticeably compressed - it has a higher peak in- 
tensity and narrower width compared to its noninteract- 
ing counterpart. This is because the internal soliton gets 
extra confinement from the ring-like soliton surrounding 
it. A similar phenomena, known as "phase separation", 
was first predicted in two-component BEG [l3, HH, in 
spherically-symmetric trap. 

The powers Nj as the functions of soliton parameter A 
at cr = —0.3 arc given in Fig. [21 It is easy to understand 
that at A = 1 two equations of the set (O, ([7]) coincide, 
that is why the diagrams A^i(A) and A^2(A) meet at A = 1. 

Let us introduce an effective radii i?„ and Ry of the 
soliton components as follows: 

r + oo r+oo 

Rl = 2ttN-^ I u^r^dr,Rl = 2ttN-^ i v^r^dr. 
Jo Jo 

As is seen from Fig. [31 the effective radius and ampli- 
tude of the u-component tend to the finite values, while 
for w-component the width increases rapidly and the am- 
plitude may come close to zero at the limit A — > 0. 
This observation explains why here the power of the u- 
component, which forms the central core, coincides with 
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FIG. 1: (Color online) Examples of stationary vector soliton solution at fixed soliton parameter (A = 0.5) and different values 
of coupling constant a. Shown are the radial profiles u(r) (solid black line) and v{r) (dashed red line). 




FIG. 2: (Color online) Beam powers (solid black line) 
and Nv (dashed red line) vs soliton parameter A for coupling 
constant a — —0.3. The insets give the examples of radial 
profiles u{r) (solid black line) and v{r) (dashed red line) for 
the points indicated on the diagrams Nu{X) and Nv{X). 



the threshold power of single 2D fundamental soliton |28| : 
Nu Nth ~ 11.7 at A ^ 0. Indeed, as A approaches 
zero, the envelope w-component practically vanishes at 
the whole localization region of u-component. 

The amplitudes (solid curves) as the functions of A 
are compared with the values of the radial profile at the 
beam axis (dashed curves) in Fig. [31 The splitting of 
the curves U = max(u) and u(0) indicates that the ra- 
dial profile u(r) gets the local minimum at r = 0, thus it 
deviates from gaussian-type shape for A > 1. The same 
effect of spatial separation of soliton components due to 
the extrusion of w-component occurs for A < 1 (see also 
the insets in Fig. [5]). This symmetry follows directly 
from the definition of the soliton parameter A — /3i//32- 
We have observed that increasing of repulsive interac- 
tions leads to the steep shrinking (within narrow limits 
in the vicinity of A = 1) of the region where both com- 
ponents are bell-shaped. Furthermore, there is no sym- 
metric vector solutions at A = 1, if cr < — 1. It is easy 
to see that for this case the set (|6]), ([7]) degenerates into 
single NLS equation that have no localized solution at 
a < —I. However, at A 7^ 1 there are non-symmetric 
"phase-separated" steady-states. This peculiarity gives 
rise to a bifurcation at A = 1 in the Nj{X) diagrams, if 
CT < — 1, as in the example in Fig. 2] (b). 
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FIG. 3: (Color online) Numerically found stationary states 
at (T = —0.3 (a) Amplitudes U = max(u), V = max(i;) and 
values at soliton center r = vs soliton parameter A. (b) 
Effective widthes of the soliton components vs A. 



B. Variational analysis 

The results of the numerical simulation can be illus- 
trated through the variational analysis. A common vari- 
ational procedure, which usually gives a good analytical 
approximation for stationary solutions, fails to describe 
a state of vector solitons with spatially separated com- 
ponents. Indeed a fixed profile of the trial function is 
not able to catch the strong modification of the spatial 
distribution observed in numerical solutions of station- 
ary NLS equations. Moreover, there is no vector solitons 
with symmetric Gaussian-type profiles in both compo- 
nents for cr < — 1. Thus, an appropriate variational pro- 
cedure should include a possibility for modification of a 
soliton shape. 



4 



We introduce a trial function with variable radial pro- 
file of the form: 

^,(r)=A,{l + <5,(r/a,)2}e"^^, (8) 

where Sj > are additional variational parameters that 
describe the modifications of the soliton profile. It is easy 
to sec that the j-th component gets a local minimum at 
r = if (5j > 1/2. The amplitudes Aj can be excluded 
using the normalization conditions 1^. Thus we have 
four variational parameters: aj and Sj {j — 1,2), which 
are found from the condition that the stationary solution 
corresponds to the extremum of the Hamiltonian at fixed 
beam powers Nj . The results of variational analysis are 
given in Figs. |4] at cr = —1.1. As it should be, one and 
only one of the parameters Si or 62 is not equal to zero for 
each approximate solution since if one soliton has a hat- 
like intensity distribution, then the other component has 
the maximum at the center. As is seen from Fig. |4] (b) 
the results of variational analysis are in good quantitative 
agreement with our numerical simulations of the vector 
solitons. 



IV. STABILITY ANALYSIS 

Before reporting our findings on stability of the 2D 
vector solitons, we review briefiy the previous results on 
this subject. The sufficient condition for collapse of mul- 
ticomponent vector solitons is found in Ref. as fol- 
lows: H < 0, where H is the Hamiltonian. This rule 
follows from the generalized on multicomponent systems 
well-known virial relation (see e.g. 0]). For the station- 
ary states H = Q, which means that the localized wave 
packet close to the stationary state can be collapsing so- 
lution. Indeed, in focusing Kerr media the 2D vector soli- 
tons are linearly unstable, as was demonstrated in Ref. 
[2^ where the generalized stability criteria similar to the 
Vakhitov-Kolokolov criteria @ is derived. Obviously, the 
repulsive intercomponent interaction has a stabilizing ef- 
fect, thus whether this type of vector solitons is stable or 
not remains to be seen. 

We solved numerically the dynamical equations (P) and 
([2]) initialized with our computed vector solutions with 
added gaussian noise. Numerical integration was per- 
formed on the rectangular Cartesian grid by means of 
standard split-step fourier technique. The snapshots of 
typical unstable evolution of the two-component bright 
soliton are given in Figs. [5][71 In these figures the in- 
tensity distributions in (x, y) plane of both components 
at different z are shown in grayscale: the darker region 
corresponds to higher amplitudes. 

Clearly, a simultaneous collapse of the both soliton 
components is not possible. To illustrate this we consider 
evolution of the stationary solution with soliton param- 
eter A < 1 when both components have close values of 
power: Ny > Nu > Nth- The snapshots of the intensity 
distributions l^'jP in ix,y) plane, found by numerical 



simulation of dynamical equations with initial conditions 
^'i(r, 0) = u{r), 'I'2(r, 0) = v{r), are given in Fig. [5] 
for A = 0.95 and a = —0.3. From the start we have the 
two-component soliton with bell-shaped intensities l^'j p. 
The 'I'2-component (which in the present case has greater 
power than the other component) gradually shrinks to 
the beam axis while the ^'i component gets the hole in 
the intensity distribution because of soliton-soliton repul- 
sive force. 

The structures exhibiting the most promise for sta- 
ble propagation arc wave packets when over-critical TV > 
Nth ring-shaped component traps the bright lower-power 
collapse- free component with N < Nth ■ The examples of 
such solutions are given in Fig [2] (the insets A and C). 
In fact, the ring-soliton keeps from spreading the inter- 
nal bright soliton in an effective potential well. At the 
same time, the repulsive core is expected to arrest the 
collapse of ring-shaped component. In our simulations, 
we indeed observed essential stabilization of such vec- 
tor solitons. However, the internal component gradually 
leaks out of the potential trap. This ever so slow tunnel- 
ing of the internal field is followed by smoothing of the 
effective potential well. As the result the bright high- 
power beam with maximum at the center is appearing 
instead of the initial ring-shaped soliton. Such an over- 
critical wave packet, obviously, is unstable with respect 
to collapse. The example of described evolution is shown 
in Fig. [51 Though we observed the robust propagation 
over hundreds of diffractive lengths, the vector solitons 
are not completely stable in this regime. 

For further stabilization the intensity and width of the 
effective potential well should be increased to prevent 
the internal component from spreading. Unfortunately, 
the new restriction on the way to complete stabilization 
of the vector solitons appears. The point is that if the 
power of ring-shaped component exceeds some critical 
value Ncr, the symmetry-breaking modulational instabil- 
ity develops. The initial ring decay into two filaments 
which drift off the center and collapse, as in the example 
in Fig. [71 To estimate the critical power Ncr, a simple 
rule can be used. First, we note that from the momen- 
tum conservation follows, that the number of filaments 
is not less than two. It was observed previously for the 
vortex solitons (26j , that each filament that appears dur- 
ing modulational instability has a power above the value 
necessary to create a single 2D fundamental soliton Nth, 
which yields the following rule: Ncr ?i '^Nth- This rough 
estimate is found to be surprisingly good approximation 
for the critical power of modulational instability observed 
in our numerical simulations. The modulational instabil- 
ity was previously found for scalar higher-order solutions 
such as vortex solitons [s^l , solitons with nodes [3l| and 
for soliton- vortex vector pairs [13, H^] ■ Of special note is 
fact that we reveal here the modulational instability for 
the ground state radially-symmetric solution. 
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FIG. 4: (Color online) (a) Variational parameters and 5j 
vs soliton parameter A at a = —1.1 (b) Beam powers Nu{X) 
(solid black line for variational results, black line with open 
circles for numerical results) and Nv{\) (dashed red line for 
variational results, red line with filled squares for numerical 
results) at cr = —1.1 



10 

5 



-5 

-10 

10 

5 





fPJ' z = 
• 


z = 22.5 

f 


z = 37.5 

o 






• 



-10 -5 5 10 



-10 -5 5 10 -10 -5 5 10 



FIG. 6: Evolution of intensities |>I'ip (upper row) and |'I'2p 
(lower row) of vector soliton for A — 0.1, a — —0.25. 
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FIG. 7: Evolution of intensities {"^if (upper row) and |'I'2 
(lower row) of vector soliton for A — 0.5, a — —2. 
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FIG. 5: Evolution of intensities I^I/ip (upper row) and |'I'2 
(lower row) of vector soliton for A — 0.95, a = —0.3. 



V. CONCLUSIONS 

In conclusion, the general properties and stability of 
2D soliton-soliton vector pairs in Kerr nonlinear media 
with focusing internal and defocusing cross-interactions 
arc studied. Stationary solutions are investigated by 
means of numerical modelling and approximate varia- 
tional method. It is found that for strong repulsive in- 
teraction between the two solitons their shapes change 
substantially. The field of stronger beam is squized out 
and forms the ring-like shell, while the weaker component 
is noticeably compressed - it has a higher peak intensity 
and narrower width compared to its noninteracting coun- 
terpart. 

We have undertaken extensive numerical simulations of 
(2-|-l)D dynamical set of NLS equations to study stabil- 
ity of the obtained stationary vector solitons. The differ- 
ent scenarios of instability have been observed depending 
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on the beam power of each component. If the power of 
one sohton component exceeds the double power of the 
fundamental Towns soliton: A^i > Ncr ~ 2Nth, such a 
ring-shaped beam exhibits the azimutlial modulational 
instability. As the result it decays into two collapsing 
spikes. If both soliton components have the powers be- 
low critical value Ncr, the modulational instability does 
not develop. In this case more powerful beam collapses 
and extrudes the field of the weaker beam outward from 
the center due to the repulsive cross-interaction. 

Finally, the quasi-stable regime occurs if the beam 
powers of the vector soliton components satisfy the con- 
ditions: Ni < Nth < N2 < 2Nth- Though the trapped 
weak beam gradually leak out through the potential well, 
which is formed by the envelope soliton, this process is 
very slow. The collapse of the over-critical beam does not 



develop until a significant portion of the trapped energy 
washes away and the repulsive core smoothes. 

On the one hand, our results impose a serious restric- 
tion on the power of a robust soliton-soliton pair, and 
on the other hand, they offer the new prospects for the 
experimental observation of long-lived two-dimensional 
vector solitons in two-component BEG and in plasmas. 
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